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ABSTRACT 

We study the linear perturbations of collisionless near-Keplerian discs. Such systems 
are models for debris discs around stars and the stellar discs surrounding supermassive 
black holes at the centres of galaxies. Using a finite-element method, we solve the 
linearized collisionless Boltzmann equation and Poisson's equation for a wide range of 
disc masses and rms orbital eccentricities to obtain the eigenfrequencies and shapes 
of normal modes. We find that these discs can support large-scale 'slow' modes, in 
which the frequency is proportional to the disc mass. Slow modes are present for 
arbitrarily small disc mass so long as the self-gravity of the disc is the dominant 
source of apsidal precession. We find that slow modes are of two general types: parent 
modes and hybrid child modes, the latter arising from resonant interactions between 
parent modes and singular van Kampen modes. The most prominent slow modes have 
azimuthal wavenumbers m = 1 and m = 2. We illustrate how slow modes in debris 
discs are excited during a fly-by of a neighbouring star. Many of the non-axisymmetric 
features seen in debris discs (clumps, eccentricity, spiral waves) that are commonly 
attributed to planets could instead arise from slow modes; the two hypotheses can be 
distinguished by long-term measurements of the pattern speed of the features. 

Key words: stellar dynamics, methods: numerical, galaxies: kinematics and dynam- 
ics, galaxies: nuclei, planets and satellites: formation, protoplanetary discs 



1 INTRODUCTION 

Debris discs are planetesimal discs that are detected through 
thermal infrared emission or scattered starlight from dust 
formed in recent planetesimal collisions. The bolometric lu- 
minosity from detectable debris discs is typically > 10~ 5 of 
the stellar luminosity, the inferred dust masses are typically 
< 1M©, and the ages of the host stars range from 10 Myr 
to 10 Gyr (see |Wyatt|2008| for a review). 

A variety of features in debris discs have been inter- 
preted as evidence for planets. These include structures in 
the ft Pictoris disc, including a warp (iHeap et al. 20001, 



a system of tilted rings (Wahhaj ct al. 2003) and a bright 



clump (|Telesco et al.|2005 1 ; clumps in the discs around Vega 



( Wyatt|2003[), e Eridani ( |Greaves et al.|2005[), r? Corvi jWy- 



att et al. 20051, and HD 107146 (Corder et al. 20091; the 



eccentricity of the discs around HR 4796A and Fomalhaut 



(Telesco et al. 2000 Kalas et al. 20051; spiral structure in 



the disc around HD 141569 (Clampin et al.|2003l; and sharp 



inner or outer edges in the discs around Fomalhaut and HD 
92945 ( |Kalas et al.|2005| |Golimowski et aL|2011[ ). 

Detailed dynamical models have shown that most or all 
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of these features can be produced by planets (see Wyatt 2009 
for a review). Moreover, in the case of j3 Pictoris (Lagrange 
et al. 20101, and perhaps Fomalhaut (Kalas et al. 20081, 



planets have been detected that may indeed be responsible 
for some or all of these features. Nevertheless, it is important 
to ask what long-lived structures could arise in debris discs 
without planets. 

In this paper we examine the possibility that low-mass 
discs can support long-lived normal modes maintained by 
the self-gravity of the disc. Normally it is assumed that de- 
bris discs cannot support such modes because of their small 
masses; all localized disturbances are dispersed by the Kep- 
lerian shear. However, a special feature of Keplerian orbits is 
that eccentric orbits do not precess. Thus the evolution of ec- 
centric disturbances in a debris disc is governed by the non- 
Keplerian forces, however small these may be. In this paper 
we shall focus on the non-Keplerian forces arising from the 
self-gravity of the disc. We neglect other possible perturba- 
tions for a variety of reasons. We ignore gravitational forces 
from planets because our principal goal is to understand the 
properties of discs in the simplest case, when no planets are 
present. We ignore radiation pressure, even though this af- 
fects the dynamics of the dust that dominates the thermal 
infrared emission and the scattered light; our justification 
is that the large planetesimals that generate the dust are 
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unaffected by radiation pressure but we recognize that the 
distribution of (invisible) parent bodies and (visible) dust 
is likely to be different. We ignore gas drag since old debris 
discs contain little or no gas, and since the planetesimals are 
likely to be large enough to be insensitive to drag. We ignore 
collisions between planetesimals because they are likely to 
be rare; indeed such collisions probably drive the long-term 
erosion of the disc in which case the collision time cannot 
be much less than the disc age. 

Debris discs are distinct from protoplanetary discs: the 
latter are comprised mostly of gas, not dust or planetesimals; 
they are much younger (typically less than a few Myr) and 



more massive (0.001-0.1 Mq) than debris discs (see Williams 
& Cieza 2011 for a review). Protoplanetary discs are de- 
pleted by various processes, including photoevaporation, ac- 
cretion onto the host star, condensation of refractory ele- 
ments into dust grains and then planetesimals, and stellar 
winds. Eventually they are likely to evolve into planetesi- 
mal/debris discs. Although our analysis here is restricted to 
collisionless systems, many of our results — in particular the 
existence of stable, slow, lopsided modes supported by the 
self-gravity of the disc — also apply to protoplanetary gas 
discs and may explain some non-axisymmetric features of 
these discs. 

To summarize we treat debris discs as collisionless sys- 
tems composed of particles influenced only by the gravity of 
the host star and the self-gravity of the disc. Their dynamics 
is therefore similar to the dynamics of discs of stars orbit- 
ing near the supermassive black holes found in the centres 
of most galaxies. Examples of these include the disc(s) of 
young stars found in the central parsec of the Milky Way 
(Genzel et al. 2010[ ), the two discs — one of young stars at 
~ 0.1 pc and one of old stars at ~ 1 pc — found at the centre 
of M31 (Bender et al. 20051, and the stellar discs that are 



inferred to form in the outer parts of quasar accretion discs 
( Goodman|2003 1 

The properties of the normal modes of low-mass near- 
Keplerian discs were investigated by Tremaine (2001; here- 
after T01), who found that (i) the frequency of the mode 
is proportional to the ratio p of the masses of the disc and 
central star, but the shape of the mode is independent of /i 
so long as < 1 (hence these are called 'slow' modes); (ii) 
all slow modes are stable; (iii) in discs with rms eccentricity 
e rms <IC 1 all slow modes have azimuthal wavenumber m = 1, 
i.e., they are lopsided. 

The results in T01 are based on linear normal-mode cal- 
culations for discs composed of particles in circular orbits, 
with softened self-gravity used to mimic the effects of the ve- 
locity dispersion or non-zero eccentricities of the particles. 
These calculations are supplemented by analytic results us- 
ing the WKB (short-wavelength) approximation, which as- 
sumes that the wavelengths of the normal modes are small 
compared to the radius. The WKB results appear to provide 
a useful guide even though this short-wavelength approxima- 
tion is not realistic for some of the disc modes. In this paper 
the effects of the velocity dispersion are computed directly, 
and we examine discs with a range of rms eccentricities e rms , 
from nearly zero ('cold' discs) to ~ 0.35 ('warm' discs). Our 
numerical results are derived using a finite-element method 
(FEM) for studying the linear normal modes of collisionless 



self-gravitating discs, as described in Jalali (2010). In partic- 
ular, we intend to address the following questions: (i) What 



are the properties of the frequency spectra of near-Keplerian 
discs? (ii) Are there any unstable modes? (iii) Are there iso- 
lated oscillatory modes in the spectrum that survive Landau 
damping? (iii) What are the differences between the spectra 
of cold and warm discs? (iv) How can stable density waves 
be excited in such discs? 

We introduce a family of axisymmetric near-Keplerian 
discs in i|2]and construct their equilibrium phase-space dis- 
tribution functions (DFs) in ij3] We obtain the governing 
equations of the perturbed dynamics in ^4] and explain the 
numerical solution procedure in ij5] We present the fre- 
quency spectra of our discs in ^ and discuss the charac- 
teristics of eigenmodes in warm and cold discs. We describe 
how these waves can be excited by tidal forces in EjT] The 
reader who is mainly interested in the application of our re- 
sults to debris discs and galactic nuclei can focus on Figures 
[7| and [TO] and the discussion in Sj8] 



2 THE MODEL 

We introduce a simple model of annular discs around mas- 



sive objects by subtracting two Toomre (1963) discs with 



n = 1 and n = 2; the resulting surface density is 
3M d f 1 1 



S d (r) = 



4tt& 2 (Jl + (r/&) 2 ] 3 / 2 [1 + (r/6) 2 ] 5 /2 
3M d (r/b) 2 

4tT& 2 [1 + ( r /fe)2]5/2' 



(1) 



where M<j is the disc mass, & is a length scale, and r is the 
radial distance to the central star. The potential correspond- 
ing to the surface density Sd is 



*d(r) 



GMd 1 + 2 (r/b) 2 
2b [1 + (r/6) 2 ]3/2' 



(2) 



with G being the gravitation constant. For a central star of 
mass M*, the total potential governing the motion of parti- 
cles is 

GM+ 



+ $d(r). 



$o(r) = 
We define 

fl= l£' R = r/b > 
and work with the dimensionless unperturbed potential 

6$o 1 U 1 + 2R 2 



Vo(R) 



GhU 
and density 

b 2 S d _ 3a* 



R 2(l + i? 2 ) 3 / 2 



R 2 



(3) 



(4) 



(5) 



(6) 



E °(' R )- Mi, 47T (1 + i? 2 )5/2 

The top panel in Figure [I] shows the radial profile of £o/M- 
The velocity of circular orbits, v c (R), is determined 
from 

2 (m _ R dV _1 u R 2 (2R 2 ~1) 
v c (R)-R— -- + - {1 + R2)5/2 ■ (7) 

The second term on the right side of |7f becomes negative 
for R 2 < 1/2. This means that our discs cannot exist in the 
absence of a central point mass. More precisely, v 2 > at all 
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Figure 1. Top: Surface-density profile of the composite Toomrc 
disc (eq. Bottom: Variation of the precession rate Q pI for p = 
0.1. For radial orbits (eccentricity e = 1) Q pT = 0. 



R if and only if p, < 5 5 ^ 2 ; this is not a limitation in practice 
since protoplanetary discs are expected to have (i< 1. 

We restrict ourselves to razor-thin discs since the verti- 
cal structure of thin discs should not strongly affect their 
large-scale response. Using the polar coordinates (R, (f>) 
and their corresponding generalized momenta {pR,p<f,), the 
Hamiltonian function governing the motion of particles 
reads 



H {p R , P4n R) = E = 



P% 



+ i>+V (R). 



(8) 



Since is a cyclic coordinate, its conjugate momentum P( f, 
is a constant of motion in the unperturbed disc. The orbital 
energy E is another integral of motion. Canonical perturba- 
tion theories describing the motion of particles, and the per- 
turbed collisionless Boltzmann equation, are substantially 
simplified by using the action variables J — (Jr, J^) and 
their conjugate angles w — (wr,!!),/,) with 



Jr = — j> PR dR, J = ~ <f> Pcj, d(f> = P4> . 



(9) 



These integrals are taken along the orbits, which consist of 
slowly precessing Kepler ellipses when p, <C 1. The unper- 
turbed Hamiltonian Ho depends only on the actions, not the 
angles. The action J^, = L is the magnitude of the angular- 
momentum vector. In the angle-action space, the equations 
of motion become 



, dHo(J) j 
w=Sl(J) = dj > 3 



0. 



(10) 



and the orbital frequencies J~2(J) = (Qr, are computed 



from 
2tt 



dR 



fi^(J) _ Jrf, 



dR 



fifl(J) I P r(r,J)' n R {J) 2tt J r 2 Pr {r,j)' 

In the limit p, — > 0, the potential is Keplerian and we have 
Q.r — = a _3//2 with a being the orbital semi-major axis. 
For < /i< 1 the radial and azimuthal frequencies are no 
longer equal, but their difference f2 pr = fi^ — £Ir is small, and 
that is the precession rate of the line of apsides. The Taylor 
expansion of Q pl begins with terms of C(/x) because fi pr 
vanishes for Keplerian orbits. Consequently, for /j« 1, the 
precession rate is proportional to the disc mass. For nearly 
circular orbits, the precession rate is given analytically by 



3/x 7? 3/2 (l-4i? 2 



(12) 



4 (l + J? 2 ) 7 /2 

Instead of the actions one may use the semi-major axis 



a (J) and eccentricity e(J) denned by 

^min (J) + -R max (J) _ -Rmax(J) 



-Rmin(J) 



2 Rmin (J) + -R max (J) ' v ' 

where 7? m in( J) and i? ma x( J) are the minimum and maximum 
distances of particles from the central star. These definitions 
are consistent with the standard Keplerian definitions when 
the disc mass vanishes. In the bottom panel of Figure [l] 
we have plotted the variation of Q pl versus a for fi = 0.1 
and several choices of e. It is seen that the precession rate 
of orbits — of any eccentricity — has a positive peak within 
the region where Eo is rising, and then switches sign and 
remains negative in the outer regions. The precession rate 
crosses through zero near a — 0.5 at all eccentricities. The 
maximum precession rate for nearly circular orbits and n <C 
1 is given by equation (12 l as uiq — 0.05861/i, which occurs 



at R — 0.2859. In !|6j we shall show that the pattern speeds 
of stable waves are closely related to uiq. 



3 PHASE-SPACE DISTRIBUTION FUNCTION 

Particle orbits in collisionless discs are not necessarily cir- 
cular. We therefore construct phase-space distribution func- 
tions (DFs) that enable us to distribute non-circular orbits 
in the disc. We seek DFs of the form ( Sawamura||1988 Pi- 
|chon fc Lynden-Bell|1996l ) 



}o{S,L) = L 2K+2 g K {£), £ = -E, 



(14) 



where < L < L c (£), L c (£) is the angular momentum of a 
circular orbit with energy E — —£, and K is a positive in- 
teger. To reproduce the surface density the DF must satisfy 
the relation 



MR) 



d£ 



fo(S,L) dL 
£)Y'\ 



L 2 



-V , 



(15) 



where L max = i?[2(* - £)] L,Z . Substituting fl4] into |15l 
and performing the integral over L we find 



B K+- 



MR) 



(16) 



2' 2) J Q R 2K+2> 

where B(p,q) is the beta function. Taking the (K + 2)th- 



order derivative of both sides of ( 16 1 with respect to 9 gives 
an explicit analytic form of <?a', 



lK+2 



MR) 



■ s /Tr2 K+1 r(K + 3/2) d*^+ 2 R 2K + 2 



(17) 
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One needs to know explicitly the function R{9) before 
doing the derivatives on the right side of (16 1. Since n is 



small in the discs we are considering, we utilize a pertur- 
bation method to compute R in terms of Let us define 
u = 1/7? and rewrite (JsJ in the form 

1 u(2 + u 2 ) 



* = u + nQ(u), Q{u) = 



2 (1 + u 2 ) 3 / 2 ' 



(18) 



We now assume a formal series expansion for u in terms of 
fj, as ( |Bellman|19"64| ) 



(19) 



and substitute this into ( 18 1. The functions Uj (^) are recur- 
sively determined by putting equal to zero the coefficients 
of jj, 3 (j — 0, 1, 2, ■ • ■ ). The recursion begins with tto = 
Up to the third-order terms, we obtain 



Hi = — Q(uo), 

U2= — ItlQ'(uo), 

uz = — u 2 Q'(u ) — ^u\Q"(uo), 



(20) 



where Q'(u) = dQ/du. The series for u converges rapidly so 
keeping the terms of 0(fi 2 ) is quite sufficient for computing 
R(iff) = in discs with n < 0.1. 

The functions Qk{£) admit negative values for K — 
0,1, but they are positive-definite and therefore physical for 
plausible values of fi < 1 when K > 2. We have plotted the 
contours of log 10 (/o//i) using (a, e) as independent variables 
in Figure [2] for K = 5 and A" = 29. The mean and rms 
eccentricity of the disc, e and e rms , are given by 



Jef (J) d 2 J = 
- ffo(J) dU - 

/e 2 /o(J) d 2 J 
" J/o(J)d 2 J 



r(f)r(f +#) 



r(3 + 

2 

2K + 5 



K) 



(21) 



Larger values of K correspond to colder discs. For /i <C 1 
the mean eccentricity e = 0.329 for K — 5 and e = 0.159 
for K — 29. When K 3> 1 the DF at a given energy or 
semi-major axis approaches the Schwarzschild or Rayleigh 
DF, 

fo(e 2 )de 2 oc exp(-e 2 /eo)de 2 , 2 = K + 1/2, (22) 

In this limit the mean and rms eccentricity are related to eo 
by e = y/neo/2, e rms = e . 

A necessary condition for stability to small-scale ax- 
isymmetric disturbances is that Toomre's Q > 1; here 
Q = a_Rn_R/(3.36Eo) where ctk is the radial velocity disper- 
sion. The models in this paper with (i<l have Q > 0.5/ fj, 
everywhere and thus are stable in this sense. The top two 
panels of Figure [3] show the rms eccentricity and or as func- 
tions of radius; for (i < 1 these are independent of /i. The 
bottom panel shows \iQ which is also independent of /i for 
H <g; 1 . Note in particular that the rms eccentricity is almost 
independent of radius. 



4 PERTURBED DYNAMICS 

We assume that gas drag, collisions, and other non- 
gravitational effects are negligible so the disc can be treated 




log 10 (a) 




log 10 (a) 

Figure 2. Contours of log 10 [/o(a, e)//i] for fj, = 0.1. The max- 
imum surface density So(-R) occurs at R = 0.8165. There- 
fore, the highest phase-space density appears in the vicinity of 
log 10 (<x) ~ 0. Top: K = 5. Bottom: K = 29. 




Figure 3. The rms eccentricity, radial velocity dispersion, and /iQ 
(the disc/star mass ratio times Toomre's stability parameter Q) 
as functions of radius. When /i <C 1 all three plots are independent 
of /i; the curves are from numerical models with fj, = 0.1. 
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as a collisionless fluid. We impose small-amplitude distur- 
bances to the surface density, potential, and DF: 

E{R,<t>,t)=T, a (R)+eSi(R,<f>,t), (23) 

V(R, 4>, t) =V (R) + eVi (R, 0, t) + eV e {R, (f>, t), (24) 

f(w,J,t)=MJ) + ef 1 (w,J,t), (25) 

where e < 1 and V c is an external perturbing potential, 
perhaps induced by a binary companion, an encounter with 
a passing star, or the tidal field of the birth cluster. The 
perturbed surface density Ei and its corresponding potential 
Vi are related through Poisson's integral: 

Zi(R',<l>',t)R' dR' d(f>' 



Vi(R,<j>,t) 



-G 



If 



y/R 2 + R' 2 - 2RR' cos((j} - ft) 
+ QR J J Si (Rf, </>',t) cos(0 - <p') dR d<j>' ^ (26) 

and we consider self-consistent density perturbations so that 



/i d v. 



(27) 



The second term on the right side of ( 26 1 is the indirect 



potential perturbation that arises because we are working in 
a non-inertial reference frame centred on the star. It is non- 
zero only for m — 1 perturbations since perturbations with 
m / 1 leave the centre of mass of the disc unchanged. For a 
particle with actions J, the radial distance R and exp(im</>) 
can be expanded as Fourier series in the angle variables, 



(28) 



Any function of R and <j> that is 27r-periodic in <f> can thus 
be expressed in the (w, J) coordinates. For the Hamiltonian 
function that governs the motion of particles, we write 

H = "Ho (J) + eVi(w, J, t) + eV e (w, J, t). 

where Ho is defined in equation |l 
turbed equations of motion become 



(29) 

Therefore, the per- 



FYU r) 



dJ 



J = 



dU 
dw 



dJ 

= -e^-(V 1 +V e ). 
ow 



(30) 
(31) 



It is obvious that the actions vary slowly in the perturbed 
disc. Subtracting the evolutionary equations of w$ and wr 
gives the apsidal precession rate in the perturbed disc, 

( d d \ 
W4,-W R = Q pl (J) + e — - 1Pr (Vi + v c ) . 



(32) 



Since Q pr = O(fi), for low-mass discs (fi 1) this equation 
contains two small parameters, /i and e. 

The DF in the perturbed disc obeys the collisionless 
Boltzmann equation (CBE), 

% = % + VM=°, (33) 

where [•, •] denotes a Poisson bracket. Here we confine our- 
selves to the linearized equation: 



dt 



1 r + [fi,n ] + [f ,v 1 ] = -[f ,v e \. 



(34) 



The remainder of this paper is devoted to the study of so- 
lutions of this equation and their application to collisionless 
near-Keplerian discs. 



5 THE FINITE-ELEMENT METHOD 

The dynamics and stability of collisionless discs are usually 
studied by one of two numerical methods: (i) iV-body simu- 
lations (e.g.,|Sellwood 1987]); (ii) expansion of the perturbed 
gravitational potential in a set of basis functions, followed 
by the evaluation of a matrix representing the response of 
the disc to a given imposed potential (e.g., Kalnajs 1977). 
Neither of these methods, however, is ideal for investigation 
of the oscillations and response of low-mass near-Keplerian 
discs, for several reasons: (i) slow oscillations are stable 
(T01) and therefore more difficult to detect than growing 
modes; (ii) slow oscillations have low frequencies, and thus 
iV-body simulations must be followed for many dynamical 
times; (iii) low-mass discs also support short-wavelength fast 
(i.e., frequency independent of /i) oscillations and these can- 
not be resolved without a large set of basis functions; (iv) 
we shall find that some slow oscillations have nearly singular 
components. Here, we adopt a finite-element method (FEM) 
and reduce the linearized CBE to a system of ordinary dif- 
ferential equations that describes the temporal evolution of 
the disc, both the eigenfrequency spectrum of an isolated 
disc and the response of a disc to external perturbations. 
We use a Co FEM (all functions are continuous, but not 
necessarily differentiable at boundaries between elements) 
in the configuration space. 

In this section, we briefly review the principles of FEM 
modelling. For a general introduction see Zicnkiewicz et 



al. (20051. Detailed descriptions of the application of an 



FEM to collisionless self-gravitating systems can be found in 



Jalali 



( 2010 1 for perturbed systems and in Jalali & Tremaine 



(2010) for equilibrium models. 



5.1 Finite ring elements in the configuration space 

We split the configuration space into TV ring elements. The 
nth element is characterized by its nodes at R n and R n +i, 
and by a linear interpolating vector G n (R) defined by 



G n — Gl 



G-> 



G hn = 1-R, G 2 ,n = R, (35) 



where R = (R — R n )/ AR n and AR n = R n +i — R n - Since 
we are interested only in linear perturbations, disturbances 
of different azimuthal wavenumber m are independent. For 
the wavenumber m, the potential V\ and the surface density 
Ei are thus computed from 



Vi (R, 0, t) =Re J2 H n (R) Gn (R) ■ a n (t)e [r 

n=l 
N 

Ei (R, <t>, t) =Re J2 H n (R) Gn (R) ■ b n (t)e ir ' 



(36) 



(37) 



The function H n (R) is unity for R n < R < R„+i and zero 
otherwise. The column vectors 

On = [ a„i a„2 ] , b n = [ b„i b n 2 ] , 

contain the nodal potentials and densities, respectively. 
According to the definition of G„(R), Ei is equal to 
Re b n i exp(jm0) at R = R n and to Re b„2 exp(im</>) at 
R = Rn+i- Similarly, the nodal potentials at these radii 
involve o„i and a„2- The perturbed surface density and its 
corresponding potential are continuous and differentiable in- 
side elements and the continuity of these functions at the 
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boundaries of elements (nodes of rings) implies 

a-n.2 = dn + 1,1, b n 2 = &n + l,l- (38) 

This means that for a given m we have N t = N + 1 inde- 
pendent nodal potentials/densities. 

The angle-action representation of the perturbed poten- 
tial Vi reads 



Vi(w,J,t)=Re Y, hi,i(J,t) e i(lWR+mw *\ 

l = — oo 

where 

JV 

hi,i(J,t) = ^2<&i{n,J)-On(t), 



(39) 



(40) 



*z(n,J) = ^ f H n (R)G n e im ^-^e- [lw n dw R . (41) 

The external disturbance 14 can also be expressed in terms 
of angle and action variables. To compute the perturbed DF 
fi(w, J, t), we use Fourier series of angle variables and write 

JV oo 

h(w,J,t) = ReJ2 E E l (n,J)-£?(t)e i « w * +mw *\ (42) 

n — 1 J — — oo 



where 

Ei(n, J) = [ E n (n,J) E l2 (n,J) 



(43) 



is an interpolating vector in the action space (to be specified 
in [ 5.3 1 and 



(44) 

is a column vector of to-be-determined DFs whose elements 
should satisfy the continuity condition 

z" 2 =4 n+1) - (45) 

Equation ( |42[ ) calculates the distribution of perturbed orbits 
based on their passage through ring elements in the config- 
uration space. If an orbit stays only inside the nth element, 
its DF becomes 

f n (w, J,t) = Re Y, E i( n > ^ ■ *f (*) e^* 4 ™^. (46) 

l= — oo 

In general, eccentric orbits may visit more than one ring 



element. The summation over n in (421 takes this behavior 
into account. 



5.2 Projected evolutionary equations 



We use the conditions ( 38 \ and assemble the nodal densities 
b n (t) and potentials cin(t) in the global iVt-dimensional vec- 
tors d(t) and p(t), respectively. Similarly, zf(t) are collected 



in zi(t). We now take the inner product of (34 1 with 

e- i< - l ' WR+mw «> ) [E l ,{n',J)] T , 

and integrate the resulting systems of equations over the 
angle-action space to obtain the Galerkin-weighted residual 



form of ( 34 1 as 



Ui(0 



dt 



-iU a (0 • zt(t) + iU 8 (Q ■ p(t) + iZ^t). (47) 



Here Ui, U2 and U3 are constant square matrices of dimen- 
sion N t x N t , and Zi(t) is an Nt -dimensional column vector, 
which is the Galerkin projection of —[fo, V e ]. 

One can also verify that the Galerkin projections of ( 26 1 



and ( 27 1 respectively become 



p(t) = C-d(t), d(t)=£)F(!)-*(t). 



(48) 



The constant matrices C and F(Z) are of dimension Nt 



N t . We combine (47 1 and (481 to express p(t) in terms of 
zi(t), and transform (471 to a non- homogeneous ordinary 
differential equation for zi(t): 

dzi{t) 



dt 



= -iUr'W ■ U a (i) ■ zi(t) +iU^ x (Z) ■ Zi{t) 

+ 00 

+ J2 iUr'W-UsW-C-FCi') •«*'(*), 



(49) 



for I, I' = 0, ±1, ±2, ■ • • . By defining 

z(t) = [ ... zl 2 zl x z T +1 z\ 2 ... ] T , (50) 

and collecting the elements of (I) ■ Zi{t) (for all I = 
0, ±1, ±2, • ■ ■ ) in the global forcing vector F(t), the system 
( |49[ ) can be cast into the standard form of linear evolutionary 
equations: 

(1 



z(t) = -iA ■ z(t) + iF(t). 



(51) 



dt 

In the absence of external disturbances, F(t) = 0, the cor- 
responding homogeneous equation admits a solution of the 
form z(t) = exp(—iu)t)c that yields the linear eigensystem: 

A ■ c = uc (52) 

We find the spectrum of w using Hessenberg transformation 
of A followed by QR factorization. The eigenvector conjugate 
to a given eigenfrequency ujj is then computed using the 
singular value decomposition 



A-Ujl = Vi 1 W V 2 , 



(53) 



where W is a diagonal matrix whose elements are the sin- 
gular values of A — Wjl, and I is the identity matrix. The 
column of V2 corresponding to the smallest singular value 
is the eigenvector as"' associated with ujj . 

5.3 Interpolating functions in the action space 

In our Co FEM analysis, the local interpolating vector func- 
tions G„ (R) (also known as shape functions) can reconstruct 
the spatial profile of any oscillatory wave whose wavelength 
is sufficiently large compared to the sizes of elements. How- 
ever, we must also interpolate /1 in the action space, which 



requires defining the interpolating vectors Ei(n, J) (eq. 42 1 



To do this we use arbitrary dynamic solutions of the lin- 
earized collisionless Boltzmann equation, which should be 
an adequate representation of the DF for the purposes of 
interpolation. In the angle-action space, and using equation 
( 39 1, one can show 



e im4, G n (R) = V 1 (n,J,w)= *i( n > J) e llWR+imw * . (54) 

1= — 00 

We assume dfi/dt = —vyfi, substitute Vi(n, J, w) into the 
linearized CBE and solve the resulting equation to obtain 
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the interpolating vector in action space (cf. eq. 42 I 



, n Idfo/dJR + mdfo/dJ^ 
Mn, J) = , ri — — n — 



(n,J). 



(55) 



The physical eigenfrequency uj will thus be equal to uj c + 7 



with w c being the computed eigenvalue of ( 52 1. Varying 7 



tests the robustness of our numerical methods since the re- 
sults should be independent of 7. Our tests show that in gen- 
eral our results are insensitive to variations of 7. Neverthe- 
less, the choice 7 > max[fl pr (J)] offers better performance, 
particularly in colder models. We originally used 7 = 0, 
corresponding to the use of static CBE solutions as interpo- 
lating vectors in action space, but with this choice we found 
occasional spurious growing modes. 

We remark that we do not generate a finite-element 
mesh in action space, for two reasons: (i) to reduce the size of 
the Galerkin-weighted evolutionary equations; (ii) to avoid 
creating spurious growing modes. The second of these prop- 
erties has a straightforward mathematical explanation: the 
number of reachable eigenmodes in the configuration space is 
equal to the number of independent nodal variables, which is 
N + l in our FEM analysis. However, the eigenvalue problem 
( 52 1 has been formulated in the phase space and the number 
of computed eigenmodes is equal to (iV + l) x (Z max — /min + 1) 
where i m i n < and Z max > are the lower and upper bounds 
in the I sums. Since the nodal densities d are related to z\ 
through equation (481, there will be (N + 1) x (i max — Imin) 



computed eigenmodes more than N + l eigenmodes that the 
dimension of d determines. The extra modes should there- 
fore overlap in groups of (7 max — Imin) members to avoid 
spurious modes. This happens in our numerical calculations 
performed in S|6]when the Fourier expansions over wr con- 
verge inside all ring elements, as is expected for the recon- 
struction of Vi and /1 in the (w r , «;0)-subspace. Generating 
a finite-element mesh, let us say with N a nodes in the two- 
dimensional J-space, will result in 2N a x (Z max — Imin + 1) 
modes, but assuming the convergence of Fourier series, only 
(N + l) groups of them will correspond to eigenmodes in the 
configuration space. Consequently, 2N a — N + l computed 
modes will be spurious, and our calculations show that such 
spurious modes are growing. Working with 2N a = N + 1 
will not help because it does not necessarily guarantee the 
convergence of FEM model in the action space. 

Only few modes out of N + 1 possible states in the con- 
figuration space (see i|6| are physical. The rest are either sin- 
gular, or do not satisfy the boundary conditions as R — > 00. 
Note that for frequencies uj that lie between the maximum 
and minimum of the precession frequency £2 pr the singular 
modes may be van Kampen modes (restricted to the surface 
in action space on which uj = fl pl ) which would be damped 
by the Landau mechanism. However, not all modes with 
frequency in this range are necessarily van Kampen modes, 
since the mode may not be produced by the orbits associ- 
ated with that resonance. Thus discrete modes may overlap 
in frequency space with continuous modes. 



6 PROGRADE WAVES 

Our finite-element mesh is uniform in log radius, 

R n =io-^+^y^ N ), 



y(n,N) 



+ 



1 



2(N+1) N + l 



(56) 
(57) 



The numerical computation of the Fourier coefficients 
^fi(n, J) and then the interpolating vectors Ei(n, J) needs 
a mesh in the (a, e) space. Such a grid is not arbitrary be- 
cause there must be at least one orbit that visits the nth 
ring element in configuration space. We fulfill this require- 
ment using the following two-dimensional grid 

(04,6*) = [Ri,y{j,M e )], 

where the grid points along the a-direction exactly coincide 
with the boundary nodes of the mesh in the configuration 
space, and there are j = 1,2, •■• , M e + 1 grid points in 
the e-direction. A circular orbit is indeed assigned to each 
boundary of a ring element. This is particularly helpful in 
cold discs where one must interpolate the population of cir- 
cular orbits. The parameters ai and 02 are chosen so that 
the computed disc mass 4tt 2 J f(J) d J using the grid points 
in the (a, e)-space agrees with the actual disc mass within 
1%. 

In this paper we focus on slow modes with azimuthal 
wavenumber m — 1. Slow modes exist with larger m, so 
long as the calculation includes Fourier terms with index 
I = — m. In particular, we have found a number of isolated, 
non-singular m — 2 modes; these are present only if we use 
a fine FEM mesh, since they are more compact and have 
shorter wavelengths than the m — 1 modes. The wavelengths 
of m = 2 modes shrink to zero as the disc becomes colder 
(see Appendix). This behavior is expected since the only 
large-scale slow modes in cold low-mass discs have m = 1. 
We found no unstable modes, which is also expected for low- 
mass discs (T01). 

We began our calculations with = M e = 70 and 
I = —1, and increased the number of Fourier terms and 
ring elemen ts u ntil the eigenfrequencies of stable modes 
found from (521 converged to a fractional accuracy of 10 -4 . 



Typically this required computing all Fourier terms with 
-2 < I < 3 and a grid with N = 160 and M e = 140 
(N = 180 and M e = 140 for the models with the lowest 
rms eccentricity, corresponding to K — 29). We have also 
experimented with including terms with larger values of \l\ 
but these had only a small effect on our results. Taking grid 
points in the regions with tiny values of f(a, e) (see Figure[2| 
leads to large errors in the properties of the calculated den- 
sity waves because the FEM discretization errors become 
larger than the absolute magnitudes of physical quantities. 
We evade this difficulty by generating the FEM mesh only 
in the annular region 0.01 < R < 100 using the parameters 
(01,02) = (2,4) inequation (57 1. 

All non-singular eigenmodes with m = 1 were found to 
be prograde (uj > 0). We find two general types of modes: a 
parent family that is already present when only the I = — 1 
Fourier component is included in the calculation, and a child 
family that bifurcates from the parent family as more Z-terms 
are included. The eigenfrequencies of child modes are very 
close to those of their parent mode (typically within 1-2%). 
They emerge from resonant interactions between two ap- 
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Figure 4. Eigenfrequency spectra of stable, prograde density 
waves in near-Keplerian discs. Only parent modes are shown. The 
vertical axis is the mean eccentricity e and the horizontal axis is 
CD = Uj/uiQ where uiq = 0.05861/j + 0(fj, 2 ). These models corre- 
spond to DFs of the form {lljl with K = 5, 10, 20, 29. The calcula- 
tion includes Fourier terms I = —2, —1, 0, 1. Note the logarithmic 
scale of the horizontal axis. Top: fi = 0.025 and uq = 0.00146. 
Middle: fi = 0.05 and loo = 0.00293. The eigenfrequencies of 
modes Di and D2 are very close and indistinguishable in the 
plots. They are CDr>i = 4.659 and CDu 2 — 4.647. Similarly, we have 
(Djjj = 4.700 and CDn 2 = 4.689. Note the similarity of the spectra 
in the top and middle diagrams despite the change of a factor of 
two in the disc mass fi; this feature is characteristic of slow modes. 
Mode shapes associated with the labelled frequencies have been 
plotted in Figures [5] [6] and [7] Bottom: Eigenfrequency spectra de- 
rived from the WKB approximation described in the Appendix. 
Each plotted point represents a degenerate leading/trailing pair 
of modes. 



proximate modes that are weakly coupled: the parent modes 
and singular van Kampen modes. For I — and / = +1 
the singular components of the child modes correspond to 
the corotation (CR) and outer Lindblad (OLR) resonances, 
respectively. The coupling between slow and van Kampen 
modes is probably due mostly to highly eccentric orbits that 
are perturbed by the gravity from both waveforms. The main 
evidence for this is that as the mean eccentricity e shrinks, 
child modes collapse to singular modes and disappear. 

We denote the maximum precession rate of circular 
max[fipr(J)]; 



orbits by loo = max[f2 pr (J)]; from equation (121 loo = 
O.O586lAt+0O 2 )- We then plot our results using the normal- 
ized frequency Co = lo/loo- Figure [4] shows the eigenfrequency 
spectra of prograde m = 1 parent modes for the mass ratios 
H = 0.025 and /j, = 0.05, and for four values of the mean ec- 
centricity e. The frequencies of child modes are not shown to 



avoid overcrowding the diagrams. Although the maximum 
precession rate loo is proportional to [i, the spectra of lo agree 
to within 1% in the models with /j, — 0.025 and fj, = 0.05. 
This scaling shows that our results can be directly applied 
to all discs with mass ratios fi <C 1, in particular to the tiny 
mass ratios fi < 0(1O~ 3 ) characteristic of debris discs. 

Figure [4] shows that the modes become more closely 
spaced as their frequency decreases and the minimum fre- 
quency in each spectrum is an accumulation point. This im- 
plies the existence of prograde waves with arbitrarily short 
wavelengths. There is also a nice correlation between the 
precession rate of the most eccentric orbit in the model (see 
Figures [I] and [2J and the lowest frequency in the spectrum. 
Models with highly eccentric orbits have an accumulation 
point of lower frequency. Figure [4] shows that the number 
of modes increases with decreasing e. In the limit of e — s> 0, 
however, dispersion-supported waves (or p-modes) can not 
exist according to WKB results (T01). 

The frequency spacing between modes B\ and B2 is 
larger than the spacing between C\ and C2, which in turn 
is larger than the spacing between D\ and D2 (which is 
so small that the two points are indistinguishable in the 
Figure). Similar behavior is seen in the F, G, and H families 
in the middle panel of the Figure. In the limit e — > 0, the 
parent modes tagged with the numbers 2j + 1 and 2j + 
2 — 0, 1,2, •■■) become degenerate. In the language of 
T01, they form a degenerate leading/trailing pair of p-modes 
(see Appendix). The pairing process begins from modes with 
highest pattern speeds, for the resonant cavities of those 
modes are fed mostly by near-circular orbits, which are the 
only population used in WKB analysis. The child modes 
of degenerate pairs also disappear because their supporting 
eccentric orbits disappear as e — > 0. Modes with id — ¥ 1 and 
sufficiently large e engage highly eccentric orbits and thus 
lead to more complex dynamics. Eccentric orbits are indeed 
the backbones of discs, and when perturbed, they affect a 
vast radial domain while near-circular orbits have only a 
local influence on developing patterns. 

We now examine the shapes of the modes. After finding 
to, we calculate its corresponding eigenvector c, and use this 
to compute the nodal potentials p and nodal densities d 



from ( 48 1 . Defining 



X(R) 



Y(R) 



ReY,H n (R)G n (R)-bn, 

n=l 

N 

lm^2H n (R)G n (R)-K, 



(58) 



(59) 



one can compute the perturbed density patterns 

Ei(iJ,0,t) = X(R)cos(m(p- cot) -Y(R)sin(m<j>-uit), (60) 

for a single wavenumber m. Note that b n are extracted from 
the elements of d using the following formula: 

K = [ d„ d n+1 ] T , n=l,2,-.- ,N. (61) 

Figure [5] shows the profile of X(R) for the labelled par- 
ent modes of Figure [4] Not only are the normalized frequen- 
cies of the modes Aj and Ej identical, but also their mode 
shapes are very similar. These remarks apply to the pairs 
(Bj,Fj), (Cj,Gj) and (Dj,Hj) as well, and demonstrate 
that the waveforms are independent of fj, so long as fj, <g; 1, 
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Figure 5. Perturbed density components X(R) (cf. eq. |60| ) for some stable modes in near-Keplerian discs of ^ = 0.025 (solid lines) 
and fi = 0.05 (dotted lines). In several panels the dotted line is not visible because it lies under the solid line. There are N = 160 ring 
elements in the configuration space for K = 5, 10 and 20, and N = 180 elements for K = 29. Filled circles mark the locations of element 
nodes. In all panels, the maximum of X(Ji) has been normalized to unity. 



as one would expect for slow modes. The figure also shows 
that the wavelength of oscillations increases with the pattern 
speed U) in a given disc, and decreases as the mean eccen- 
tricity of the disc shrinks. The number of nodes increases 
as the frequency decreases. An interesting property of the 
waves showing multiple nodes is that their density peaks are 
approximately equally spaced in logarithmic scales. 

The child modes are hybrid modes that inherit the fea- 
tures of their parents in the central regions of the disc, but 
have a spike at the location of singular modes that cou- 
ple to them. Figure [5] displays the parent mode Ds of fre- 
quency lu — 2.083 (see Figure|4| and its children Ds,cr with 
Q = 2.0765 and Dg.oLR with uu — 2.0728, which contain sin- 
gular van Kampen modes at the corotation and outer Lind- 



blad resonances, respectively. In low-mass discs, these reso- 
nances are at large radii where the surface density is small, 
so the singular component of a child mode involves only a 
small fraction of the mass involved in the parent mode. As 
the disc mass shrinks to zero the child modes merge with 
the parent mode. The reason is that the eigenfrequency of 
the parent mode is proportional to the disc mass so with 
very small disc masses the corotation and outer Lindblad 
resonances are at extremely large radii where the surface 
density is negligible. Thus the distinction between parent 
and child modes is unimportant for low-mass discs such as 
debris discs. 

Figure [7] displays shaded contour plots of the pattern 
of Ei(_R, 4>, t) for some models with fj, — 0.025 (mode shapes 
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Figure 6. Perturbed density components X(R) for the par- 
ent mode Dg and two of its associated child modes Dg and 
Dg oLRi m which the parent mode is coupled to singular modes 
at the corotation (CR) and outer Lindblad (OLR) resonances. 
The model parameters are fi = 0.025 and K = 29. 



corresponding to /i = 0.05 are similar). It is seen that the 
wave packets are more radially compact in the colder (K — 
29) model than warmer (K = 5, 10) ones. 

The properties of these modes can be explored using 
the short-wavelength or WKB approximation described in 
the Appendix. The validity of this approximation requires 
k > h/R where k is the wavenumber and h is a dimensionless 
number of order unity. If two adjacent nodes of a wave are at 
Ri and R2 then k dR — ir so the condition for validity 
of the WKB approximation may be written 7r > h log R2 / Ri 
or log 10 R2/R1 < 1.36/ft. Inspection of Figure [6] shows that 
for h = 1 this condition is satisfied by all of the modes we 
have computed, though not by much in some cases. 

The bottom panel of Figure [4] shows the WKB fre- 
quency spectrum. Each point corresponds to a degenerate 
pair of modes, one composed of leading spiral waves and the 
other of trailing. These modes arise from waves in the res- 
onant cavities defined by the closed frequency contours in 
Figure |A1| The WKB approximation correctly reproduces 
several striking features of the FEM frequency spectra: (i) 
all modes are prograde (u > 0); (ii) both the number of 
modes and the maximum frequency grow as the mean ec- 
centricity e of the disc shrinks; (iii) there is an accumulation 
point of modes near cj/ojq = 1 in the FEM spectra and at 
cj/ojo = 1 in the WKB spectra; (iv) there is also reasonable 
quantitative agreement between the frequencies derived by 
the two methods, at least for the discs with the lowest mean 
eccentricity. The WKB analysis in the Appendix fails to find 
the child modes, for two reasons: (i) it is based on the epicy- 
cle approximation, which assumes that the eccentricity is 
small and thus neglects the highly eccentric orbits that cou- 
ple the slow and van Kampen modes; (ii) it is based on the 
approximation that the disc mass /i — > 0, and in this limit 
the pattern speed of the slow mode goes to zero so the outer 



Lindblad and corotation resonances are at very large radii 
where the disc surface density is negligible. 



7 EXCITATION OF OSCILLATORY WAVES 

Protostars live in the harsh environments of their birth clus- 
ters. Simulations of the Orion nebula ( Scally fc Clarke|2001 \ 
show that about 10 per cent of stars can have encounters 
closer than 100 AU within 10 7 years. Such encounters can 
excite waves in planetesimal/debris discs. Encounters were 
invoked as a possible explanation for the asymmetries in the 
f3 Pictoris debris disc by |Kalas fc Je witt ( 1995J and |Larwood| 
& Kalas (2001 1 but these authors treated the debris disc as 



a collection of test particles, which can give misleading re- 
sults since the self-gravity of the disc dominates the apsidal 
precession. 

Since our goal here is only to illustrate this process we 
confine ourselves to in-plane parabolic encounters. Consider 
a disc particle orbiting around a star of mass M*, and as- 
sume a perturber of mass M p . As in earlier sections, we scale 
all lengths so that the disc length scale b is unity, and denote 
the normalized position vectors of the particle and perturber 
(with respect to the host star) by R and R p , respectively. 
The equation of motion for a disc particle is 



= -V[a k -R + V (R) + Vi(R,t) + V P (R, R p )] , (62) 



d 2 R 

where a* is the acceleration vector of the host star in 
an inertial frame, Vo is the unperturbed gravitational po- 
tential due to Mi, and the self-gravity of the disc, Vi is 
the perturbed self-gravitational potential of the disc, and 
V p = -(M p /M i ,)/\R P - R\ is the potential field of the per- 
turber. The gradient V is taken over the R space, and the 
normalized time t is related to the actual time i ac tuai through 

Victual = (GAU/b 3 ) 1 ' 2 . 

We assume that a* is due to the encounter; thus we 
ignore the cluster's tidal field. Consequently, 



R 



Mp \ R v R 



R% 



Rp — \Rp\ 



(63) 



For a distant encounter, R <C R p , the potential V p can be 
expanded as the following series 



Vo 



Mp 
M t 



1 (<p - <f> p ) 



1 

R p 
R„ 



E 

i=0 

R 



Pi [cos (0 - (p p )] . 



(64) 



R P R 



where <fi and <f) p are, respectively, the azimuths of the disc 
particle and perturber measured from an inertial reference 
line, and Pi are Legendre polynomials. The effective poten- 
tial due to flying-by perturber thus reads 



V c = a^-R + Vp, 

'Mp\ 1 ^ 
M-, R v ^ 



Pi [cos ( 



(65) 



where we have dropped the i — term in ([64]) because it 
makes no contribution to the force, and the i = 1 term has 
been cancelled by a+ ■ R. 

Modes having azimuthal wavenumber m = 1 can only 
be excited by those i > 3 terms of V e that produce cos <j> 
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Figure 7. The patterns of oscillatory waves in the configuration space for a near-Keplerian disc with the mass ratio fj, = 0.025. We have 
displayed only the positive part of £i(B, <j>,t) at t = 0. Maximum densities of all panels have been normalized to unity, and the contour 
levels range from to 1. The point mass sits at (0, 0). Left: K = 5. Middle: K = 10. Right: K = 29. 



and sm(j> factors. Modes with m = 2 are excited by the 
i = 2 term of V p , which is much larger than the i > 3 terms 
for distant perturbers (R/R p <C 1); however, we have found 
(*j6| that slow modes with m = 2 have wavelengths that are 
generally smaller than those of m = 1 modes, even for discs 
with a relatively large mean eccentricity e, and which shrink 
to zero as e — > 0. Thus m = 2 modes couple less effectively 
to smooth perturbing potentials. We conclude that the dom- 
inant slow mode excited by an external perturber may have 
either m = 1 or m = 2. 

For brevity, we shall examine only m = 1 modes here. 



The dominant term of ( 65 1 for m — 1 perturbations is 

y. - - - 1 i ( — i ( .: I cos(0-<^). (66) 



This can be expressed in the angle-action variables as (cf. 
eq.p| 



V e ~Re Q(t)hc,i(J)e K 



where 



Q{t) 



1 



i0 p (t) 



(67) 



(68) 



MJ [R P (t)V 

is the time-varying part of the external perturbation, and 
1 



h c ,i(J) = 



2tt 



R cos [Iwr + (wj, — (f>)] dwR. 



(69) 



For a parabolic encounter with minimum distance i? p , m i n 
and gravity parameter M — 1+Mp/M*, the true anomaly <j> p 
and radial distance R p are computed through the following 
equations: 



t{<t> p ) 



R« 



V2 



l.u. I ^ 



+ ^tan 3 ^ 
3 V 2 



2R„ 



1 + cos( 



M 



2 R 3 



(70) 
(71) 



The effect of the external perturbation on the evolution 
of /i inside the nth element is determined by the Galerkin 



projection of the Poisson bracket — [/o, V e ] as 
47T 2 iZl5'(t) = - J! E?,(n',J)[f ,V B ] 



(72) 



Substituting from ( 67 1 into ( 72 \ and performing the integral 



over the angle space gives the 2-vector 



ZT(t) = Q(i) 



8Jr dJ<b 



h e> i(J)Ef(n,J) d 2 J, (73) 



whose components (Z" l and ; ) are, respectively, the con- 
tribution of the disturbing force to the inner and outer nodes 
of the nth ring element in the configuration space. The dis- 
turbance at the jth ring node thus reads 



^1,1 > 



ryN 

\ ^2,1 J 

and we obtain 
Zi(t) = [ Zx,i Z 2 ,i 



3 = 1, 

Kj<N + l, 
j = N + l, 

Zn,i Z(jv+i),i 



(74) 



(75) 



Defining Fi = Uj x (i) • Zi, the global forcing vector is assem- 
bled as 



F(t) = [ ... i^(t) F^(t) F T +1 (t) 
= Q(t)g, 



(76) 



where g is a constant vector. For each Fourier number I, 
we have N t unknown DFs collected in the vector zi(t). The 
Fourier series in terms of wr is usually truncated at some 
i m in < and i max > 0. We thus have M = (Z max — imin + 
1) x TVt unknown DFs that we collect in the N- vector z(t). 
Similarly, F(t) and g are A/"-dimensional vectors. 

Any excited wave is a superposition of all eigenmodcs 
of d52l: 



z(*)=5>(i)z W , 
J=l 



(77) 
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Figure 8. The forcing vector components Cj for the modes of the 
discs having fi = 0.05 and K = 20 and 29. The corresponding 
frequency spectra have been plotted in Figure [4] 



but not all eigenvectors z^' are physical. Increasing the 
number of elements increases the accuracy of the eigenvalues 
and eigenvectors describing the isolated oscillatory modes 
but also adds spurious and/or singular modes. Such non- 
physical modes can contribute noise to the calculated disc 
response. To keep only physical modes, we introduce 



31 12 



m = r 



„(A0 



and express ( 77 \ in the matricial form z = M • q. This is 



q = iJ ■ q+iQ^M" 1 



(78) 



substituted into (51 1 to obtain 

d_ 
dt 

where J = M 1 ■ A - M is a diagonal matrix — or a Jordan form 
if there are degenerate eigenvalues ( Perko|2001 1 — whose el- 
ements are the eigenfrequencies of (521. The diagonalizing 



matrix M is often called the modal matrix. 

We call the A/"-dimensional vector £ = M~ ■ g the forc- 



ing vector and rewrite ( 78 1 in terms of its components: 



^Ij (*) = -iwjgj (t) + iQ(t)Q , j 



1,2, 



(79) 



Equations associated with non-physical cjj can now be 



dropped from ( 79 \ and we find both the homogeneous and 



particular solutions, 

qj(t) = e-^ (t - to) gj -(to) + iO f Q(r)e- iw ' ( '- T) dr, (80) 

for j = 1,2,... , A/" p with M p being the number of physical 
modes. Fly-by perturbations begin at t = to = — oo {cj) p = 
— 7r) with qj(to) = (Vj). Consequently, using the orbit 
equations ( |70[ ) and ( |71[ ), and defining fi = u p /ujj, we arrive 
at 



3i(*>0p) = - 



3^2(M-1) 5/3 



64M 4 / 3 
0P (l + cosO 2 ^^*^^ d£, 



(81) 
(82) 



which leaves behind the permanent oscillation 

... .3V2(M-1) 5/3, n , m -tayt 

when the encounter ends at p = +7r. The integrands of the 
real and imaginary parts of Qj(n,/3) are, respectively, even 



(83) 
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Figure 9. The profile of Qj (w, /3) (solid line) together with its 
analytic asymptotes when /3 — > (dashed line; eq. |85[ l and /3 — > oo 
(dotted line; eq. |84[ l. The excitation is efficient for /3 > /3 C with 
/3 C = 0.169367 (see text for definition). 



and odd functions of 4> p over the interval [— 7r, +tt]. One thus 
obtains Im[Qj(7T,/3)] = 0. The real part of Qj(n,(3)/l3 5/3 is 
positive-definite, and therefore, a necessary and sufficient 
condition for the excitation of the jih oscillatory mode is 
that the corresponding component of the forcing vector Q 7^ 
0. The asymptotic forms are 

0,(^/3) «2^/3 5/3 , 



for p > 1 and 



q15/4„1/2 
2 ' 7T ' 1/6 

p ' exp 



2%/2 \ 
" 3/3 J 



(84) 



(85) 



for p <C 1. The exponential decay for small w p is due to 
adiabatic invariance. 

We have computed £ for all parent modes of Figure 
[4] and have plotted its components Q in Figure [8] for two 
H = 0.05 models with different mean eccentricities. The re- 
sults for child modes and other models are similar. In our 
models Q is larger for modes with low frequencies. This can 
be understood as a competition between two effects seen in 
Figure [5] (i) as the mode frequency decreases, the number of 
its nodes increases, so the coupling of the mode to a smooth 
external field is reduced; (ii) as the frequency decreases, 
the outermost peak of X(R) shifts to a larger radius and 
hence contributes more to the term — [/o, V e ] in (34 1 — recall 
that V c ~ R 3 . In general the second effect wins, causing the 
coupling, as measured by £, to be larger for low-frequency 
modes. 

We also find that the range of £ is similar for all K 
models (Fig. |8J. This shows that the response of a near- 
Keplerian disc is not sensitive to its mean eccentricity: m = 1 
slow modes in warm and cold discs have an equal chance of 
being excited by encounters. 

The excitation efficiency of modes is determined by the 
function Qj(7T,P), which has been plotted in Figure [9] The 
excitation of mode j is inefficient for uj p < p c ujj where we 
have defined the critical frequency ratio /3 C — 0.169367 as the 
point where Qj(n, 0) drops to 10% of its value predicted by 
its fi — > 00 asymptote (eq. 841. For a given mass parameter 
M, a perturber can only excite mode j efficiently if its orbit 
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Figure 10. The evolution of mode Di in a disc with mean eccentricity e = 0.159 as a perturber on a parabolic orbit encounters the 
star-disc system. <j> p is the azimuthal angle of the flying-by star; the line 4> p = coincides with the x-axis and with the direction of 
pcriastron. We have <p p = (— 7r,+7r) for t = (— oo,+oo). The contours show the positive part of the response density. The periastron 
distance of the perturber is /? p , min = 188.58(1 + M p /M*) 1 / 3 (M d /10~ 3 A/*)~ 2 / 3 . 
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has periastron i? p , m i n < [M / {P^uf)] 1 ' 3 . Faster modes have 
larger ujj , and therefore need closer encounters to be excited. 

In Figure [To} we illustrate the excitation of mode Di at 
four azimuths during the fly-by of a disc with mass ratio fi — 
0.025, for an encounter with the parameter lj p — \/2^Di = 
0.0096538. The interaction begins at <f) p — —n (t = — oo) 
and ends at <f) p — +tt (t — +oo). We plot the positive part 
of the response density: 

E 1>Dl (i?,<M) = Re{gj[t,tfp(t)] e 1 "*-"*** 

N 

xY,H n (R)G n (R)- (86) 

n = l 

where j corresponds to mode Di and the vectors &„ are 
extracted from d^') = £\ F(0 ' z i as we did in equation 

i el). 

We remark that stable modes always rotate with a con- 
stant angular velocity, but the perturber-star centreline has 
a variable angular velocity. Therefore the perturber may lead 
or lag the maximum response in azimuth, and the maximum 
response may occur some time after closest approach. 

In general, of course, the close passage of a per- 
turber will excite multiple modes. The main visual 
difference between a single-mode and multi-mode re- 
sponse is the occurrence of long-period beating pat- 
terns in the latter case. We have constructed anima- 
tions of the evolution of the multi-mode pattern dur- 
ing an encounter and the beating can be quite strik- 
ing ( |http://www.youtube.com/watch?v=ZTyXK7H 6Q8E) 
This animation is for a model with K — 10 and jj, — 0.025, 
and the 11 modes with the highest frequencies are partici- 
pating in the response. 



8 APPLICATION TO DEBRIS DISCS AND 
GALACTIC NUCLEI 

We have shown that low-mass, near-Keplerian, collisionless 
discs can support stable, long-lived, large-scale slow modes. 
The most prominent of these are expected to have azimuthal 
wavenumbers m = 1 and m — 2. 



8.1 Debris discs 

The existence of slow modes implies that debris discs can 
support waves in the planetesimal population that pro- 
vides most of the disc mass, and suggests that collisions 
in this non-axisymmetric distribution could generate non- 
axisymmetric dust distributions that would be visible in 
thermal emission or scattered light. 

Non-axisymmetric structures in debris discs are nor- 
mally assumed to be produced by planets, but our results 
imply that some or even most of these structures may be 
density waves. Specific examples include: 

• /3 Pictoris: Scattered starlight reveals that this star is 
surrounded by a debris disc extending to > 1000 AU. The 
disc is brighter on one side than the other, perhaps due 
to an m = 1 slow mode, and also contains brightness en- 
hancements that could be due to shorter- wavelength density 
waves. The disc exhibits warps or tilted rings at various radii 
( Heap et al.|2000||Wahhaj et al.|2003[ ); although the present 



paper examines only in-plane slow modes, there should also 
be slow bending modes, and these provide a possible expla- 
nation for the warps. There is a 10Mj planet orbiting at 
~ 10 AU in the f3 Pic system (Lagrange et al.|2010 1 but it is 
far from clear that this is the cause of the warps and other 
features; several authors have argued that the asymmetries 



provide evidence for two or even three planets (Frcistetter 
et al.|2007 Currie et al.|2011 \ but it is implausible to invoke 



a new planet for every feature 

• Fomalhaut: This star is surrounded by a ring of dust 
with a sharp inner edge at 130 AU. The centre of the ring is 
offset by 15 AU from the host star, implying an eccentricity 
of 0.11; the ring is narrowest at apastron, implying that the 
eccentricity declines with radius ( Kalas et al.|2005" I . Quillen 
(2006) stressed that these features could be produced by a 
planet orbiting just inside the ring, and a possible planet was 



subsequently discovered (Kalas et al. 20081. The eccentricity 



of the ring could be forced by the planet or a slow density 
wave, depending on whether the planet mass or ring mass 
is larger. The sharp inner edge of the ring is most likely due 
to the planet. 

• Vega: Observations at a variety of wavelengths between 
350 /im and 1.3 mm reveal a face-on dust ring of radius ~ 100 



AU, dominated by two clumps (see Marsh et al.||2006] for a 
summary of the data, and beware that |Pietu et al.|2011| ques- 
tion the reality of non-axisymmetric structure in the disc). 
The clumps are usually ascribed to dust trapped in a reso- 
nance with an unseen planet (e.g., Kuchner fc Holman|2003] 
Wyatt 2003]), but m = 1 and m = 2 slow density waves pro- 
vide an alternative explanation. Within a few years we may 
be able to distinguish these hypotheses by measurements of 
the motion of these clumps relative to the host star: the ex- 
pected angular speed of the planet is ~ l°yr _1 while slow 
modes should have negligible pattern speeds. 

• e Eridani: A nearly face-on ring of dust surrounds this 
star at ~ 60 AU. The disc exhibits several clumps and a 
lopsided brightness distribution in images at 450 (jm and 



850 fim ( Greaves et al. 2005 1. Some but not all of these peaks 



may be background sources. Models in which the clumps are 
due to resonances with a planet are describ ed by|Ozernoy et 
all (120001), iQuillen & Thorndike|(|2002|), andlDeller & Maddi- 



(2005j). These features could be due to slow modes, but 
the presence of density maxima at several azimuths would 
require that more than one mode was present. Resonance 
models predict angular velocities around the host star of 
about 1° yr _1 . 

• HR 4796A: There is an edge-on debris ring ~ 80 AU 
from the host star. One ansa of the ring is brighter, hotter, 



and at smaller radius than the other (Telesco et al. 2000 



Moerchen et al. 20111. This asymmetry is most naturally 
explained by an eccentric dust ring (e ~ 0.06); at perias- 
tron the dust is closer to the star and therefore hotter and 
brighter ("pericentre glow" 



Wyatt et al. 19991. The ring 



eccentricity is usually assumed to be excited by secular per- 
turbations from a nearby planet but an m = 1 slow mode 
of the disc is an alternative. The mode might be excited 
by the companion star HR 4796B, currently at a projected 
separation of ~ 500 AU. 

• AB Aurigae: Near-infrared images reveal a debris disc of 
over 1000 AU radius. The disc shows spiral arms at radii of 
several hundred AU, some of which are also seen at submm 
wavelengths, as well as rings, gaps, and clumps at smaller 
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radii (Hashimoto et al. 20111. As in the case of other sys- 



tems, the features could be due to planets or slow modes, 
and these hypotheses can be distinguished by proper-motion 
measurements. 

• r/ Corvi: The disc surrounding this star appears at 
450 /im as two equally bright peaks equidistant from the 
host star at 100 AU; these can be modelled either as the 
ansae of an edge-on axisymmetric ring or as a more face-on 
disc containing dust trapped in a resonance with a planet 
( |Wyatt et al.| [2005). A third possibility is an m = 2 slow 
mode in a face-on disc. 
141569A: 



HD 



Clampin et al. (20031 have detected 



strong spiral structure in the debris disc around this star, 
at about 400 AU radius. They suggest that the spiral may 
be excited by tides from nearby stars. There is also a gap 
in the disc at about 250 AU radius; both of these features 
might be due to planets but only if planets can form at radii 



exceeding 200 AU. Wyatt (20051 has suggested that the spi- 



ral could be caused by a Jupiter-mass planet on an eccentric 
orbit (e ~ 0.2, a ~ 250 AU) but slow modes provide a more 
economical explanation, especially given the difficulties of 
forming planets at such large distances. 

• HD 100546: This disc exhibits an apparent dark hole 



and bright clump at about 30 AU from the host star ( Quanz 
|et al.|[2011| ). These features could be due to an orbiting 
planet or a slow density wave. The Keplerian motion at this 
radius is about 3°yr _1 . At much larger radii, ~ 250 AU, 
the disc exhibits spiral structure ( |Grady et al.| [2001). Pos- 
sible explanations include a planet at several hundred AU 
from the star or density waves excited by a passing star. 
The latter possibility was discussed by Quil len et al.| ( |2005[ ) 
but dismissed because their estimated lifetime for the spi- 
ral structure was only ~ 10 4 yr and no suitable nearby star 
could be found; the results of the present paper imply that 
the structure could last for a much longer time — perhaps 
as long as the 10 Myr age of the star — so the chance of a 
suitable encounter in the past is much larger. 

• HD 61005: This star is surrounded by an asymmetric 
edge-on debris disc of radius ~ 60 AU. The asymmetry can 
be modelled as a mean eccentricity of 0.05, but there are no 
planets more massive than ~ 3 Jupiter masses close to the 
ring ( |Buenzli et al.|2010 l. 

• HD 15115: This star hosts an edge-on debris disc; the 
dominant thermal emission from the disc arises at radii ~ 35 
AU but the disc is visible to much larger radii. The sur- 
face brightness of the east side of the disc is about 1 mag 
fainter than the west side at a given radius and the surface- 
brightness distribution perpendicular to the disc midplane 
is asymmetric on the west side (Kalas et al. 20071; both 



features can arise naturally from an m = 1 distortion. 

• HD 107146: There is a dust ring at 100 AU that exhibits 
clumps and a lopsided brightness distribution in 1.3 mm 



images (Corder et al. 20091. These might be due either to 



a planetary resonance or to slow density waves; however, 
880 /im observations with similar resolution do not confirm 
the existence of the clumps ( Hughes et al.poTT I . 



8.2 Discs in galactic nuclei 

The results of this paper also illuminate our understanding 
of stellar discs in galactic nuclei. They can be applied di- 
rectly to such discs if the apsidal precession is dominated by 



the self-gravity of the disc, rather than relativistic effects or 
the gravitational field from a spherical stellar population in 
the nucleus. 

The apparent 'double' nucleus of M31 is most likely 
a stellar disc that has been distorted by a large-amplitude 
slow mode (see |Peiris fc Tremaine||2003| |Salow fc Statler| 



2004 



and references therein). Such modes arise naturally 
in N-body simulations ( Jacobs fc Sellwood|200l" l. They can 
be excited by gas inflow and star formation in the central 



few parsecs of the galaxy (Hopkins & Quataert 2010 1 or 



by instabilities induced by a small population of counter- 
rotating stars ( |Touma|[2002[ ). Slow modes may also play a 
central role in feeding supermassive black holes (Hopkins & 
Quataert|20"Tot . 



9 DISCUSSION 

The finite element formulation has enabled us to explore 
the modal spectrum of low-mass near-Keplerian collisionless 
discs, and to calculate the corresponding mode shapes for a 
wide range of initial radial dispersions (rms eccentricities). 
Our method also yields moments of the distribution func- 
tion, which provide the evolutionary equations for energy 
and angular-momentum transport in perturbed discs, and 
allows the accurate representation of modes that contain a 
singular resonant component (e.g., Fig. [6J. 

We find that near-Keplerian discs support 'slow' modes, 
that is, modes for which the eigenfrequency or pattern speed 
is proportional to the disc mass. WKB analysis shows that 
these modes are closely related to the p-modes found by 
Tremainc (2001) in cold near-Keplerian discs with softened 
gravity. Both our numerical results and analytic arguments 
imply that there are no unstable slow modes. All slow modes 
in the discs we have examined are prograde (positive pat- 
tern speed). Slow modes can exist with arbitrary azimuthal 
wavenumber m, but modes with m = 1 and m — 2 have 
the largest scale and are the easiest to excite by an external 
perturber. 

The eigenmodes of the linearized CBE bifurcate from 
the degenerate leading/trailing modes predicted by WKB 
theory. The modes are degenerate for cold discs and split 
into close (in frequency) pairs as the mean eccentricity 
grows, until for e > 0.2 there is no apparent pairing in fre- 
quency space (see Fig. EJ. 

Some of the non-axisymmetric structure that is com- 
monly observed in debris discs, such as clumps, lopsided 
rings, and spiral arms may be due to slow modes, perhaps 
excited by the fly-by of a passing star or binary companion. 
These features are normally ascribed to hypothetical mas- 
sive planets embedded in the disc. The two hypotheses can 
be distinguished in some discs by monitoring the motion of 
these features over decade timescales: many features asso- 
ciated with planets should orbit the host star at a pattern 
speed that is not far from the Keplerian angular speed of 
the planet, whereas slow modes should have negligible pat- 
tern speeds. Structures induced by modes and planets may 
also be distinguishable in the future by high-resolution far- 
infrared observations by interferometers or large single-dish 
telescopes (e.g., ALMA or CCAT). 

Future theoretical work should include the exploration 
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of slow bending modes and of the behavior of slow modes in 
thick discs that resemble the discs seen in galactic nuclei. 
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APPENDIX A: WKB ANALYSIS 

The dispersion relation for a collisionless disc can be com- 
puted analytically in the WKB or short- wavelength approx- 
imation ( Binney fc Tremaine|2008 I 

1 - mij^ / a R k \ 2 



_27rGE |fc| 



where 



G(s,x) = -e" 
X 



E T 



— s 2 /n 2 



(Al) 



Here Qr and fi^ are the orbital frequencies (eq. Ill, I n (x) is 
a modified Bessel function, m is the azimuthal wavenumber, 
(j_r is the radial velocity dispersion, and the perturbed sur- 
face density is assumed to vary as Ei (R, (f>, t ) oc e xp\i(m</> + 
f k(r')dr' — cut)]. The dispersion relation (All is valid if 
or « Q^R and \k\R > 1. 
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Figure Al. Contours of constant uj in the WKB approximation, from equation | |A2| for m = 1 (left) and m = 2 (right). The disc has 
mean eccentricity e = 0.159. The horizontal axis is the wavenumber and the vertical axis is the log of the radius (base 10). Lighter shades 
correspond to higher values of ui. 



For low-mass discs, /i < 1, the precession frequency is 



D.r = 0(jj) (eq. 12 1 and the surface density 



Eo is also O(fi). For slow modes ui is 0([i), and since s = 
(cj — mQ)/Qn is nearly — m, the denominator 1 — s 2 /n 2 
becomes small for n = m and the solution of the dispersion 
relation for a given m is obtained by keep ing the dominant 
n — m term of the summation in (All. 



s R 2 fl 2 R where e rms 



Moreover a 2 R 



is the rms eccentricity defined in 
equation ( |21[ ) — strictly, this is the rms eccentricity at a given 
radius rather than of the whole disc but in our models the 
rms eccentricity is almost independent of radius (Fig. j3j). 
The dispersion relation for slow modes then simplifies to 



m7rEo|fc| 
Or 



T m \{e ima kRy/2] 



where F m (x) = ~e x I m (x)- 
X 



(A2) 



Here all quantities are written in the dimensionless units of 
^ Note that as z -> 0, I m (z) -t z m /(2 m m\). Thus as the 
rms eccentricity shrinks to zero the wavenumber of a slow 
mode must vary as k ~ ems - "'"' 2 '" -1 '. In other words, for 
m = 1 slow modes have \kR\ ~ 1 even for cold discs — in 
this case the use of the WKB approximation for slow modes 
is not formally justified, but the results provide a useful 
qualitative guide to the behavior of the frequency spectra 
that we find using FEM (see further discussion at the end of 
§61. For m > 1 slow modes exist but with wavelengths that 
shrink as e rms declines. For cold discs J- m (x) — for m > 1 
so the disc supports only singular modes at the resonances 
UJ = rafipr. 

In the WKB approximation, disturbances in the disc 
can be decomposed into wavepackets that propagate at the 
group velocity dw/dfc along contours of constant frequency 
cj. These contours are illustrated in Figure |A1| for a disc 
with mean eccentricity e = 0.159 and m = 1,2. Wavepack- 
ets that propagate along open contours eventually wind up 
(|fe| — > oo) and disappear. Discrete normal modes can arise 



for closed contours if the appropriate resonance condition is 
satisfied. Consider the case m = 1 (left panel of Fig. All. 
For the closed contours centered on |fe| = 10, R = 0.6 the 
resonance condition is (T01 Q 

dfcdi? = 2n(n - ~) , (A3) 

where n = 1, 2, 3, . . . and the integral is taken over the area 
in (k, R) space enclosed by the contour. These modes, called 
p-modes by T01, occur in degenerate pairs, one composed of 
leading and one of trailing waves (k < and k > respec- 
tively) . 

Equation ( |A3[ ) can be solved numerically to find the 
frequencies of the p-modes. These are plotted in the bottom 
panel of Figure [4] for m = 1, and are in good qualitative 
agreement with the frequencies calculated by FEM. There 



is similar agreement between equation ( A3 1 and FEM for 
m — 2 modes. 



1 The factor of 2 on the left side of equation (56) in T01 is incor- 
rect. 



© 0000 RAS, MNRAS 000, 000-000 



